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Magnet or ot at ional turbulence in stratified shearing boxes with 
perfect gas equation of state and finite thermal diffusivity 

G. Bodo 1 , F. Cattaneo 2 , A. Mi gnone , P. Rossi 
ABSTRACT 

We present a numerical study of turbulence and dynamo action in stratified 
shearing boxes with zero mean magnetic flux. We assume that the fluid obeys 
the perfect gas law and has finite (constant) thermal diffusivity The calculations 
begin from an isothermal state spanning three scale heights above and below the 
mid-plane. After a long transient the layers settle to a stationary state in which 
thermal losses out of the boundaries are balanced by dissipative heating. We 
identify two regimes. A conductive regime in which the heat is transported mostly 
by conduction and the density decreases with height. In the limit of large thermal 
diffusivity this regime resembles the more familiar isothermal case. Another-the 
convective regime- observed at smaller values of the thermal diffusivity, in which 
the layer becomes unstable to overturning motions, the heat is carried mostly by 
advection and the density becomes nearly constant throughout the layer. In this 
latter const ant- density regime we observe evidence for large-scale dynamo action 
leading to a substantial increase in transport efficiency relative to the conductive 
cases. 

Subject headings: accretion disc - MRI - MHD - dynamos - turbulence 

1. Introduction 

The origin of turbulence and enhanced angular momentum transport in accretion flows 
is a fascinating problem of considerable importance in astrophysics. It is commonly believed 
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that the magneto- rotational instability (MRI), in one form or another, plays a fundamental 
role in destabilizing the basic quasi-Keplerian flow. When a net magnetic flux is present the 
MRI sets in a s a classical linear instab ility with a well defined growth rate and characteristic 
wavenumber ( IBalbus fc Hawleylll99ll ); the turbulence then develops from the nonlinear evo- 
lution of this instability. When there is no net flux the problem is more complicated and the 
turbulence must develop from a nonlinear subcritical instability. In this case, the problem 
becomes fundamentally one of establishing what form of dynamo action can be sustained in 
a disc. Much of what is currently known about dynamo action in accretion flows is based 
on numerical studie s formulated within the framework of the shearing-box approximation 
( IHawley et al.lll995l ). The simplest set up, both conceptually and numerically, consists of an 
unstratified, isothermal shearing-box with periodic boundary conditions in the vertical direc- 
tion. It is now well established that this configuration suffers from the so-called convergence 
problem. As the magnetic diffusivity decreases, or equivalentl y the resolution increases, the 
Maxwell stresses decrease, eventually to become negligible (IFromang fc Papaloizoul 120071 . 
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see however Fro- 

mang 2010 for a different view). The cause for this "non convergence" has been attributed to 
the the lack of a characte ristic outer scale in the periodic, unstratified problem (for a discus- 



Bodo et al.l 120 111 ). The next step towards more realistic simulations is to retain the 



sion see 

shearing-box geometry but with the inclusion of vertical gravity, and consequently, stratifi- 
cation. This introduce s a characteristic length-the scale height-that may help to remedy the 
convergence problem (IDavis et al.l |2010| . IShi et al.l |2010| . lOishi fc Mac Low! l201ll ) . Whether 
this is the case or not, at the moment, remains an open question. Certainly, in the strat- 
ified cases the sol utions manifest a richness both in space and time that is absent in the 
unstratified cases ( Gresselll2010l . iGuan fc Gammiell201ll . ISimon et al.ll2012l ). It is important 
to note that most of these studies adopt an isothermal equation of state; the resulting den- 
sity distribution is correspondingly close to a Gaussian with most of the mass concentrated 
near the mid-plane, and tenuous, low density regions above and below. This leads to very 
different dynamo processes operating in the mid-plane and in the overlying regions. Al- 
though an isothermal formulation is conceptually simple and easy to implement numerically, 
it neglects the possibly important process of turbulent heating by viscous and Ohmic dissi- 
pation. It can be argued that in an optically thin environment turbulent heating may not 
be important since the energy can easily escape without substantially heating the ambient 
plasma. However, this is definitely not the case in an optically thick environment. In this 
case the plasma will be heated locally and the final thermal structure will be determined by 
a balance between energy deposition and energy transport. In this case, it is possible that 
substantial departures may develop from the isothermal case that, in turn may impact the 
oper ation of the dynam o . Some of these issues h ave been addressed by Hirose and collabora- 
tors (IHirose et al.ll2006l . 120091 . iBlaes et al.l 1201 if ), who have considered radiation dominated 
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discs and have included a so phisticated treatment of the radiation field, and also by the works 



of Flaig and collaborators ( IFlaig et al.l l2010l |2012| ) whose models of proto-planetary discs 



include partial ionization, chemical networks and heat transport in the radiative conduction 
approximation. All these works indicate that turbulent heating can indeed be important. 
Here, we also address the problem of turbulent heating but in the somewhat simpler case 
of a fully ionized, pressure dominated disc. Our intention is to provide a bridge between 
the works of Hirose et al. and Flaig et al. and those based on the isothermal equation of 
state. To this end we consider a stratified shearing-box with a perfect gas equation of state 
and finite (constant) thermal diffusivity. The objective is to study how the basic state and 
the corresponding dynamo action changes as the thermal diffusivity is varied. In this work 
we deliberately keep the formulation as simple as possible in order to highlight some of the 
basic underlying physical processes. 



2. Formulation 

Our objective is to provide a simple model in which the effects of dissipative heating 
can be studied. In particular we want to assess how these processes, together with thermal 
transport lead to departures from the more familiar isothermal cases. We assume that 
the plasma is optically thick and approximate the radiative transport by a diffusion process 
which we model by a thermal conduction term in the energy equation. In the spirit of keeping 
things as simple as possible, and in order to capture more easily the general properties of the 
solutions, we make further simplifications, by neglecting the dependencies on density and 
temperature resulting from the diffusion approximation to the radiative transport equation, 
and assuming a constant thermal diffusivity. A more realistic treatment of the radiation will 
be considered in future work. 

We perform three-dimensional, numerical simulations of a perfect gas with thermal 
conduction in a shearing box with v ertical gravity. A det ailed presentation of the shearing 



box approximation can be found in lHawley et al.l (119951 ) . The Magneto-Hydro-Dynamics 



(MHD) shearing-box equations, including vertical gravity and thermal conduction can be 
written as: 

^ + V-(pv) = 0, (1) 



dv B-VB 1 /B 2 \ / „ 1 „ 9 N 
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^-Vx(vxB) = 0, (3) 

a 171 

— + V • [(E + P T )v + (v • B)B - WT] = 0, (4) 

where B, v, p and P denote, respectively, the magnetic field intensity, the velocity, the 
density and the thermal pressure; E is the total energy density, Pt is the total (thermal plus 
magnetic) pressure and k is the thermal conductivity. The local angular velocity Q = Qe z 
and the shear rate 

A = *™ (5) 
A ~ 2 dR { ) 

are assumed constant. For a Keplerian disk A = — (3/4)f2. The system is closed by the 

equation of state for a perfect gas: 

P = pT (6) 

where we have absorbed the perfect gas constant in the definition of the temperature. The 
thermal conductivity can be written as 

k = h -K P (7) 

where k is the thermal diffusivity, which, as discussed above, we assume to be constant, and 
the factor of 5/2 is appropriate for a gas with three degrees of freedom. 

We start our simulations from a state with a uniform shear flow, v = —2Axe y , and 
density and pressure distributions that satisfy vertical hydrostatic balance with constant 
temperature T . With these conditions the initial density has a Gaussian profile given by 

p = p exp(-fiV/2T ), (8) 

where po is the value of density on the equatorial plane. 

If the MRI develops to substantial amplitude, this initial state will be driven away from 
thermal equilibrium by the energy input from dissipative processes. The temperature in the 
equatorial regions will progressively increase and a thermal gradient will be established until 
a new equilibrium is reached whereby the energy input is balanced by thermal losses at the 
upper and lower boundaries. As we shall see, the new equilibrium can be quite different from 
the initial isothermal state and is determined self-consistently by the heating associated with 
the process of angular momentum transport by the MRI. We note here that, in our current 
formulation, we do not include viscous and Ohmic dissipation explicitly. The heating of the 
fluid occurs because of numerical dissipation together with a conservative formulation of the 
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total energy equation. The latter requires that whatever kinetic or magnetic energy is lost 
by dissipative processes it be re-introduced in the form of internal energy (heating). 



The computational domain covers the region L x x L y x L z , where L 2 
and L 7 = 6 if, where 



H, L y = ttH 



H 



(9) 



is the pressure scale height in the initial isothermal state. In the vertical direction the 
box is symmetric with respect to the equatorial plane z — 0, where gravity changes sign. 
Numerically, the domain is covered by a grid of 32 x 96 x 192 grid points. A magnetic field 
of the form 

° sin (ir) * (10) 



B 



is imposed initially where B corresponds to the ratio between thermal and magnetic pressure 
and has a value of 1600. Clearly, there is no net magnetic flux threading the box. In addition 
we introduce random noise in the y component of the velocity in order to destabilize the 
system. 

Following common practice, we assume periodic boundary conditions in the y direction 
and shear periodic conditions in the x direction. In the vertical direction, we assume that 
the upper and lower boundaries (z = ±3H), are impenetrable and stress free, giving v z = , 
dv x /dz = dvy/dz = 0, and also that the magnetic field is purely vertical, giving dB z /dz = 0, 
B x = By — 0. We should note that these conditions allow a net flux of magnetic helicity 
through the boundaries with, possibly, impo rtant consequences to the dynamo processes 



(IVishniac fc Choi 1200 ll . iKapyla fc Korpill201ll ). Finally, we assume that the boundaries are 



in hydrostatic balance, and that the temperature is constant and equal to T ; thus 



dz 



T3ptt H, 



T = Tn. 



:n) 



All simulations are carried out with the PLUTO code (jMignone et al.l 120071 ). with a sec- 
ond order accurate scheme, HLLD Riemann solver and an explicit treatment of thermal 
conduction. 



3. Results 

We now describe the development of MRI driven turbulence from an initially isothermal 
state. Hereinafter, and unless otherwise specified, when presenting the numerical results, we 
adopt Q^ 1 as the unit of time, H as the unit of length, and the mid-plane density in the 
initial isothermal state po, as the unit of density and, since H is our unit of length, we have 
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To = 1/2. Following the initial perturbations, a sub-critical instability sets in leading to 
the generation of magnetic fields and the development of turbulence. Dissipative processes 
heat the plasma driving the system away from the initial isothermal state. Eventually the 
system reaches a stationary state in which the heat generated by the turbulence is balanced 
by the heat lost through the upper and lower boundaries. Locally, the balance is between 
the volumetric heat production and the divergence of the heat flux that can arise both by 
thermal conduction and turbulent transport. The relative importance of these two processes 
depends on the value of the thermal diffusivity, which, here, is expressed in units of the 
product of the scale-height and the isothermal sound speed, i.e. it has the form of an 
inverse Peclet number. A typical evolution for a case with k = 2 x 10~ 2 can be followed 
in Fig. [TJ where we show the time history of the Maxwell stresses averaged over the entire 
computational domain. Clearly, there is a long adjustment phase lasting approximately 500 
time-units after which the system settles into a stationary state in which the stresses remain 
strongly fluctuating but with a well defined (time) average value. The corresponding thermal 
history can be assessed by inspection of Figs. [2] and [3j showing respectively, the horizontally 
averaged temperature, T(z) and density p(z), at several times. The asymptotic profiles in 
the stationary state (obtained by time averaging from t = 500 to the end of the simulation, 
t = 2000) are denoted by angle brackets. Clearly the increase in the Maxwell stresses are 
accompanied by the heating of the central regions leading to the establishment of a nearly 
parabolic profile in temperature. We note a corresponding dramatic change in the density 
distribution that evolves from the initial Gaussian profile to an almost constant distribution 
at later times. 

The development of a constant density state is somewhat remarkable, and deserves 
further investigation. At first sight it may appear as the result of a fortuitous choice of 
thermal diffusivity. As we shall see presently, this is not entirely the case. In the stationary 
state the average temperature and density are related by the condition of hydrostatic balance, 
which we write here in dimensional form, and for simplicity we only consider z > 0; 



depends on the relative magnitude of the two terms in the brackets on the RHS of ({12]) . The 
first term is a fixed linear function of z. The second-the temperature gradient-is negative 
since the layer is heated from within, but its magnitude depends on a balance between local 
heat production rate and local heat transport. To a first approximation, one could assume 
that the energy production rate should be independent of the thermal diffusivity k. This 
is not unreasonable, since the production rate is driven by turbulent dissipation, which in 
turn depends solely on the efficiency of the MRI. This being the case, the magnitude of 




(12) 



Clearly, whether the density decreases upwards, increases upwards, or remains constant 
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the temperature gradient could be made arbitrarily small by choosing a large value of k. 
Clearly, if the thermal diffusivity is huge, thermal conduction can easily transport all the 
generated heat along very shallow gradients. The temperature will be nearly constant, the 
density will rapidly decrease upward in accordance to f lT2"j) . and resemble the isothermal 
distribution. By contrast, if k is tiny the temperature gradients required to carry the heat 
will be huge (in absolute value), the RHS of (|T2|) will be positive and the density will 
rapidly increase with height. However, this configuration with a density inversion is strongly 
unstable to Rayleigh- Taylor type instabilities. The resulting overturning motions will both 
carry the heat more efficiently than thermal conduction, and homogenize the mass towards 
a constant density state. Thus we can conjecture the existence of a critical value of k = K CTit 
above which the transport is mostly conductive, the layers have a density decreasing with 
height and a stratification approaching that of an isothermal layer in the limit of large k 
(conductive states). For k <C n crit , the heat transport is mostly advective, and the density 
is approximately constant (convective states). 

Some of these ideas can be easily verified by considering a series of calculations with 
varying thermal diffusivity. The results are summarized in Figs. S] and where we show 
the steady state horizontally averaged temperature and density distributions for different 
values of k. As expected, the temperature gradient is always negative (z > 0), its magnitude 
increases with decreasing k as does the overall temperature of the layer. For large values 
of k the temperature distributions have an approximate parabolic profile and the density 
decreases upwards. 

For small values of k the temperature in the interior approaches a "tent" profile with a 
parabolic shape near the equator then a linear decrease over most of the domain and thin 
boundary layers at the edges. As k decreases the profiles move up retaining their shape but 
producing progressively thinner boundary layers. The corresponding density profiles confirm 
the establishment of a constant density state that becomes asymptotically independent of k. 
From Fig. we can estimate that the critical value of n for a transition from conductive to 
convective regimes, in this setup, satisfies 



Our conjecture that as k crosses its critical value the vertical heat transport changes from 
conductively dominated to advectively dominated can also be verified by considering the 
horizontally averaged conductive, and advective fluxes that can be defined, respectively, as 



2 x 10 



-2 



(13) 



5 dT 



(14) 



c = KP—^, 

2 H dz' 



and 




(15) 
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Their values in the stationary state for the two extreme values of k are shown in Fig. |6j 
The roles of the two types of flux practically reverse. For k = 0.12 the transport is entirely 
conductive and advection is negligible; for k = 4 x 1CT 3 heat conduction is negligible, except 
in the boundary layers, and all of the flux is carried by advection. It is interesting to note 
that near the equator where the advective flux is small-it is actually zero at the equator-the 
density displays a weak inversion. This is related to the absence of gravity near the equator 
to drive Rayleigh- Taylor instabilities. 



i 1 1 r 




time 



Fig. 1. — Time history of the Maxwell stresses averaged over the computational box for the 
case k = 2 x 10~ 2 . 

The existence of two regimes, conductive and convective, with strikingly different vertical 
structures is likely to lead to correspondingly different dynamo actions. A measure of these 
differences can be assessed by inspection of Fig. [9] where the domain averaged Maxwell 
stresses are shown as a function of time for different values of k. The corresponding curve 
for an isothermal case is also included for comparison. Clearly, the angular momentum 
transport efficiency increases with decreasing k and eventually saturates in the convective 
regime. It is natural to assume that once the heat transport is mostly advective further 
decreases in thermal diffusivity will not make any difference. What is remarkable is the 
difference between the convective cases and the purely isothermal one, with the latter being 
strikingly smaller. 

Further evidence for two distinct types of dynamo actions operating in the two regimes 
can be obtained by inspection of Fig. [71 This shows the horizontally and time averaged 



Fig. 2. — Temperature averaged over horizontal planes, T as a function of the vertical 
coordinate z for the case k = 2 x 1(T 2 . The different curves refer to different times, as 
indicated by the labels, and for comparison we plot also the time averaged distribution (T) 
in the steady state. 
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Fig. 3. — Density averaged over horizontal planes, p as a function of the vertical coordinate 
z for the case k = 2 x 10^ 2 . The different curves refer to different times, as indicated by 
the labels, and for comparison we plot also the time averaged distribution (p) in the steady 
state. 




Fig. 4. — Plot of (T) as a function of z. The different curves refer to different values of k, 
as shown in the legend. For comparison we plot also the isothermal case 
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Fig. 5. — Plot of (p) as a function of z. The different curves refer to different values of k, as 
shown in the legend. For comparison we plot also the isothermal case 




Fig. 6. — Plot of (F c ) and (Ft) as functions of z. The solid curves show (F c ) while the 
dashed curves show (Ft), the different colors refer to different values of k as indicated in the 
legend. 
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Fig. 7. — Plot of the average Maxwell stresses as a function of z for two cases with different 
t values of k. One case (solid line) is in the convective regime, the other (dashed-dotted line) 
is in the conductive regime. 

Maxwell stresses as a function of z for two cases with different values of k corresponding to 
the convective and conductive regimes. The curve for the conductive case follows the general 
trend of the more familiar isothermal calculations. The transport is largest in the denser 
central regions, steadily decreasing at higher values of z. This is in sharp contrast with the 
convective case in which the the stresses actually rapidly increase with distance from the 
mid-plane reaching a sharp maximum near the boundaries. In both cases, the corresponding 
Reynolds stresses are small and decrease steadily away from the mid-plane. The spatio- 
temporal behavior of the dynamo is also remarkably different in the two regimes as illustrated 
in Fig. [HI The two panels show space-time diagrams of the horizontally averaged azimuthal 
magnetic field as a function of z and time. The lower panel, corresponding to a conductive 
case, displays the characteristic patterns typical of the isothermal cases signaling the presence 
of cyclic activity with magnetic structures propagating from the mid-plane to the boundaries. 
In the upper panel there is no evidence for cyclic activity or pattern propagation. The 
magnetic structures form and vanish seemingly at random with no apparent characteristic 
time between field reversals. Furthermore there are events in which coherent structures form 
that extend over the entire layer. Interestingly, for earlier times, when the layer is still close 
to isothermal there is some evidence for pattern propagation. From these last two figures it 
is clear that both the transport efficiency and the amount of generated toroidal flux is much 
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time 



Fig. 8. — Space-time diagrams of average azimuthal field. The horizontally averaged value 
of By is plotted as a function of z and t. The upper panel corresponds to a case in the 
convective regime; the lower, to one in the conductive regime. The corresponding values of 
indicated. 
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higher in the convective regime than in the conductive one. 



A possible reason for this difference might be related to the influence of magnetic bound- 
ary conditions. It is well known that in unstratified shearing boxes the boundary conditions 
make a big difference to the operation of the dynamo. Periodic boundary conditions, as was 
mentioned in the introduction, lead to small-scale dynamo action and to the convergence 
problem. On the other hand, "vertical" boundary conditions, like the ones imposed here, 
lead to a much more efficie nt dynamo that appea rs to scale with the system size rather than 
with the dissipation scale (IKapyla fc Korpill201ll ). By contrast, the soluti ons in isothermal 



stratified shear i ng boxes are more insen sitive to the boundary conditions ( IDavis et al.ll2010 



Shi et al.ll201Cl lOishi fc Mac Lowll201l[ ). This most likely is because the boundaries are lo- 
cated in very tenuous regions characterized by low density and high Alfven speed. In the 
convective cases described here, the density is nearly constant as a function of height making 
the layer appear more "unstratified" . Partial support for this argument can be provided by 
looking at what type of dynamo is operating in the convective regime. Fig. [TUl shows the time 
history of the volume averaged value of B y (the azimuthal component) scaled in terms of the 
rms value of the fluctuations. Two things are worthy of notice: the average field changes sign, 
and its magnitude is comparable with-and actually it occasionally exceeds-th at of the fluc- 



tuati ons. This is strongly suggestive of the operation of a system-scale dynamo (ITobias et al. 



201 ll ) and should be contrasted with the corresponding isothermal case in which there is a 
different behavior depending on height and the ratio between average and fluctuations raises 
from about ten percent in the central region, where most of the transport takes place, to 
more substantial values in the upper and lower regions where the transport strongly declines. 



4. Conclusions 

Our main objective has been to study numerically the effects of dissipative heating and 
finite heat transport in determining the thermal structure of the layer and the efficiency 
of angular momentum transport in stratified shearing boxes with zero magnetic flux. In 
particular we wanted to compare with the more commonly studied isothermal case. To this 
end we have considered the simple case of a fluid obeying the perfect gas law and with finite 
(constant) thermal diffusivity. 

Our main result is to identify two distinct regimes: conductive and convective, corre- 
sponding respectively to large and small values of the thermal diffusivity. In the conductive 
regime, the heat generated by dissipation is transported through the bulk of the layer by 
thermal conduction, and the temperature and density have close to para bolic profiles. Thi s 



appears to be in agreement with the conclusions of the recent work by lUzdenskyl ( 120121 ). 
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Fig. 9. — Time histories of the volume averaged Maxwell stresses for different values of the 
thermal diffusivity k. The values of k are shown in the legend. 
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Fig. 10. — Time history of B y (the volume averaged B y ) in units of the rms value of the 
fluctuations. For this case k = 4x 10 -3 . 
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The convective regime is dramatically different. In these cases the heat is transported al- 
most entirely by overturning motions driven by Rayleigh- Taylor type instabilities. The 
density profile becomes flat, and the temperature develops a "tent" profile with thin bound- 
ary layers at the upper and lower boundaries. There is evidence that the "tent" profile 
for the temperature and the flat profile for the density are universal, in the sense that they 
depend solely on the properties of the turbulence and not on the values of the collisional pro- 
cesses. This last property in particular, may have important consequences for the dynamo 
processes. It appears that the dynamo can operate more efficiently in a layer with nearly 
constant density than in a corresponding layer with the same total mass and a Gaussian 
profile (isothermal case). This being the case, there is a interesting feedback effect. The dy- 
namo drives the MRI turbulence that heats the layer causing it to become Rayleigh- Taylor 
unstable, the overturning motions associated with the Rayleigh- Taylor instability homoge- 
nize the density allowing a more efficient operation of the dynamo, and so on until the layer 
settles to a universal, convective, self-regulated state. At the moment it is not clear whether 
the Rayleigh- Taylor driven motions contribute directly to a more efficient working of the 
dynamo, or they contribute indirectly by maintaining the more beneficial constant density 
state. Som e of the similariti e s betw een the dynamo properties observed here and those in 



the work of iKapyla Korpil (120111 ) in which stratification is absent suggest that it may be 



the constant density feature that is important. 

Finally, we remark on the natural extensions of the present model. There are two 
avenues that immediately come to mind. One is to include a more realistic treatment of 
the thermal transport. The obvious next step is to consider thermal diffusivities that have 
power law dependencies on density and temperature. We anticipate that this may have 
some impact on the stratification in the conductive regime, but hardly any in the convective 
regime in which all the thermal transport is mediated by flows anyway. The other is to 
consider more realistic boundary conditions, like for instance those appropriate to black 
body radiation. Preliminary results in this direction show that in the convective regime this 
choice leads to a change in the overall value of the temperature but not in its profile. Also, 
the constant density profile remains unchanged. These results however are preliminary and 
a more thorough study is needed. Also, the assumption of impenetrable stress-free boundary 
conditions should be replaced by more realistic conditions in which there is a thin transition 
layer across which the opacity changes dramatically and the fluid goes from being optically 
thick to optically thin. The problem is similar to that of matching a photosphere on top of 
a stellar convection zone. Numerically this is extremely challenging and will be considered 
in future works. However all these extensions are secondary to the issue of convergence. In 
a sense, if the dynamo ceases to operate efficiently at high magnetic Reynolds numbers all 
bets are off. There is some room for cautious optimism since the evidence so far, is that 
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the dynamo operating here in the convective regime is more likely to be of the system-scale 
type than of the small-scale type. Preliminary calculations with twice the resolutions indeed 
support this conjecture. However, in the end only a (very costly) convergence study will 
settle the issue. 
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